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ABSTRACT 

The main orbital signatures of the secular evolution of an isolated self-gravitating 
stellar Mestel disc are recovered using a dressed Fokker-Planck formalism in angle- 
action variables. The shot-noise-driven formation of narrow ridges of resonant orbits 
is recovered in the WKB limit of tightly wound transient spirals, for a tepid Toomre- 
stable tapered disc. The relative effect of the bulge, the halo, the disc temperature and 
the spectral properties of the shot noise are investigated in turn. For such galactic discs 
all elements seem to impact the locus and direction of the ridge. For instance, when 
the halo mass is decreased, we observe a transition between a regime of heating in the 
inner regions of the disc through the inner Lindblad resonance to a regime of radial 
migration of quasi-circular orbits via the corotation resonance in the outer part of the 
disc. The dressed secular formalism captures both the nature of collisionless systems 
(via their natural frequencies and susceptibility), and their nurture via the structure 
of the external perturbing power spectrum. Hence it provides the ideal framework in 
which to study their long term evolution. 
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1 INTRODUCTION 

Understanding the dynamical secular evolution of galactic 
discs over cosmic times is a long-standing endeavour, which 
recent renewed interest as their cosmological environment is 
now more firmly established in the context of the ACDM 
paradigm (Planck Collaboration et al. 2013). Indeed, dis¬ 
entangling the respective role of the cosmic environment 
(nurture) and the intrinsic structure of galaxies (nature) in 
explaining the observed physical and morphological distri¬ 
bution of galaxies is focussing much recent activities. Self- 
gravitating discs are cold dynamical systems, for which rota¬ 
tion represents an important reservoir of free energy. Some 
perturbations are strongly amplified, while resonances tend 
to confine and localise their dissipation: even small stimuli 
can lead to discs evolving to distinct equilibria. 

Modern N-body simulations now allow for both detailed 
modeling of intricate non-linear physical processes (such as 
multi-scale hydrodynamics, star formation, AGN feedback... 
see e.g. Dubois et al. 2014), but also well controlled ideal¬ 
ized numerical experiments (Sellwood & Athanassoula 1986; 
Earn & Sellwood 1995; Sellwood 2012). Such experiments 
are essential to understand how the orbital structure of a 
galactic disc may drive its secular evolution and simply take 


into account their self gravity. Yet the reliability of numeri¬ 
cal simulations to preserve the symplectic nature of the un¬ 
derlying physical system over hundreds of orbital times is 
potentially an issue which calls for alternative probes. 

In parallel, over the past few years it has been found 
that the formalism of angle-action variables (Born 1960; 
Goldstein 1950) and the so-called matrix method (Kalnajs 
1976) also provided means of accounting for the self-gravity 
in the secular equation driving such systems in the limit of 
large number of particles. They are in particular well suited 
to disentangling the intricate roles of resonances and identi¬ 
fying the orbital families driving their secular evolution. 

Two equations for secular diffusion have recently 
been revisited using these coordinates: i) the (possibly 
dressed) Fokker-Planck equation (Binney & Lacey 1988; 
Weinberg 1993; Pichon & Aubert 2006; Chavanis 2012b; 
Nardini et al. 2012; the companion paper, hereafter paper 
I, Fouvry et al. 2014b, submitted), when the source of (pos¬ 
sibly coloured) potential fluctuations is taken to be an 
external bath, e.g. corresponding to the cosmic environ¬ 
ment; ii) the Balescu-Lenard non-linear equation (Balescu 
1960; Lenard 1960; Weinberg 1998, 2001; Heyvaerts 2010; 
Chavanis 2012a) which accounts for self-driven orbital secu¬ 
lar diffusion induced by shot noise corresponding to the dis- 


2 Jean-Baptiste Fouvry, Christophe Pichon 


creteness of the system. These equations are fairly unique 
in providing a theoretical framework for the secular evolu¬ 
tion of stellar and dark matter dominated systems, and well 
suited to explain, complement and understand the results 
of these crafted numerical experiments. As an illustration of 
their versatility we will implement here the (simpler) Fokker- 
Planck equation in order to explain one such experiment, 
and postpone the implementation of the Balescu-Lenard 
equation to further investigation. 

Indeed, Sellwood (2012) (hereafter S12) has suggested, 
using a well controlled numerical setting, that an isolated 
stable stellar disc would secularly drift towards a state of 
marginal stability through the spontaneous generation of 
transient spiral structures. In its experiment, S12 evolves 
a set of particles, of increasing number, for hundreds of dy¬ 
namical times. These particles are distributed according to 
equation (4) below which corresponds to what should be 
a stable distribution. Nonetheless, S12 identifies the noise 
driven formation of ridges in action space, along very spe¬ 
cific resonant directions. Indeed, discs which are fully stable 
at a linear level may still on the long-term develop spiral 
structure that can grow to important amplitudes, eventu¬ 
ally transforming the disc into a barred-type spiral galaxy. 


We intend to show here how the formalism of dressed 
secular diffusion written in angle-action variables is able to 
predict the observed drifting process exhibited in S12. The 
direct analytical or numerical calculation of the modes of a 
galactic disc is a complex task, which has only been made 
for a small number of disc models (Zang 1976; Kalnajs 
1977; Vauterin & Dejonghe 1996; Pichon & Cannon 1997; 
Evans & Read 1998; Jalali & Hunter 2005). For the sake of 
analytical simplicity, following paper I, we will make use of 
the so-called WKB approximation (Liouville 1837; Toomre 
1964; Kalnajs 1965; Lin & Shu 1966), and assume that the 
initial distribution is well described by a tepid Schwarzschild 
distribution function, which will allow us to compute the 
gravitational susceptibility of our disc as a simple (scalar) 
multiplicative factor, and express the diffusion coefficient 
algebraically. This approximation provides insight into the 
location of the relevant resonances. Our strategy here is to 
defer to appendices as much of the technical detail as we can 
while still conveying an understanding of the overall theory 
in the main text (see also paper I). The main orbital sig¬ 
nature of the secular evolution of the tepid Toomre-stable 
tapered Mestel disc will be recovered in the WKB limit of 
tightly wound transient spirals. The relative effect of the 
bulge, the halo, the disc temperature and the spectral prop¬ 
erties of the shot noise will be investigated in turn. 


The paper is organized as follows. Section 2 introduces 
briefly the Mestel (1963) disc considered by S12 and presents 
its main features. Section 3 recaps the formalism of the 
secular diffusion equation in action-space, which is able to 
eapture the main observed features when considered in the 
WKB limit for a tepid disc. Our results for the Mestel disc 
are presented in section 4. Finally, we conclude in section 
5. Appendices A and B present respectively a rapid deriva¬ 
tion of the diffusion equation and its WKB limit; paper I 
presents a more systematic derivation. 


2 THE DISC MODEL 


Stellar discs are dynamical systems, which at leading order 
have reached a virialized state within an axisymmetric grav¬ 
itational field created not only by their own mass, but also 
by other constituents of the galaxy, mainly the inner bulge 
and the surrounding dark-matter halo. The disc considered 
by S12 is an infinitely thin Mestel disc (Mestel 1963), for 
which the circular speed — R /dR — Vq of the stars 
is independent of the radius. The stationary potential back¬ 
ground of such a disc and its associated surface density are 
given by 


MR) = vSiyi 
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where Vo and Ri are scale parameters. Let us intro¬ 
duce the actions of the system, {Jr,J(p) (Born 1960; 
Binney & Tremaine 2008). For an axisymmetric two di¬ 
mensional disc, the azimuthal action = Lz is the an¬ 
gular momentum, whereas the radial action is given by 
Jr = I/tt vr dR and encodes the amount of radial en- 

J -^min 

ergy of a star, where vr= [2(F —'0o(R)) —is inte¬ 
grated between the pericenter Rmin and apocenter Rmax of 
the trajectory. Here Jr = 0 corresponds to circular orbits. 
These action remain well defined through the secular evo¬ 
lution of the disc as cylindrical symmetry is preserved. For 
simplicity we use the epicyclic approximation to describe 
the behavior of the distribution function of the system in 
the action-space (Jr, Jc^). This approximation holds as long 
as the particles do not have too eccentric orbits. Because 
the Mestel disc has constant circular velocities, the link be¬ 
tween the angular momentum J^^ and the guiding radius Rg 
is straightforward and reads 


Jcf) — • 
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Within the epicyclic approximation, one can obtain the ex¬ 
pression of the intrinsie frequencies, the azimuthal frequency 
Q(J(^) and the epicyclic frequency given by 
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Two remarks should be made on these frequencies. First, 
within the epicyclic approximation, the two frequencies are 
only function of the angular momentum J^^ and do not de¬ 
pend on the radial action Jr. Moreover, one should note that 
we have n/Q. — v^, so that the Mestel disc is an intermedi¬ 
ate case between the Keplerian case for which n/Q. — 1, and 
the harmonic case for which n/Q. — 2. The ratio between 
the intrinsic frequencies is an important parameter for the 
dynamical behavior of the system, since it determines the 
location of the resonances and a constant ratio may intro¬ 
duce degeneracies. Using the epicyclic approximation, the 
distribution function considered by S12, takes the form of a 
locally isothermal-T)¥ or Schwarschild-DF, which reads 


Fo (Jr, J (j)) 
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where the intrinsic frequencies are given by equation (3), 
cTr is a constant dispersion describing the spread in radial 
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Figure 1. Surface density Et of the tapered Mestel disc. The unit 
system has been chosen so that Vq = G = Ri = 1. 


velocity of the distribution function, and the taped surface 
density Et in analogy with equation (1) is given by 


Et(Jc^) = 


Vf 


27r GJ(f) 
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(5) 


In equation (5), Tinner and Touter are tapering functions used 
to damp out the contributions from the inner and outer 
regions, which read 
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where u and ji control the sharpness of the two tapers. The 
inner tapering function at the scale Ri induces an impor¬ 
tant distribution function gradient at this scale, which is 
indeed responsible for the position of the peak of diffusion 
and reflects the presence of a bulge. The outer tapering 
function reflects the finite size of the disc. For the numeri¬ 
cal simulations, we used the same constants as in S12. We 
place ourselves in a unit system such that: Vq = G = Ri = 1. 
The other numerical factors are given by ar — 0.284, z/ = 4, 
/i = 5, i?o = ll-5. The shape of the damped surface density 
Et is shown on figure 1. The tapering functions have for ef¬ 
fect to reduce and turn off the self-gravity of the disc in the 
inner and outer regions. As such, it provides a fairly gen¬ 
eral class of models for more realistic discs. One can now 
look at the initial level contours of the distribution function 
represented on figure 2. One of the consequences of such a 
scale-invariant disc is that its local Toomre Parameter Q 
(Toomre 1964) is almost independent of the radius for the 
intermediate regions. Indeed, one defines Q as 




(Jr K/(^Jcj)^ 


(7) 


where in order to reduce the susceptibility of the disc, we 
suppose that only a fraction, , of the disc is self-gravitating, 
with 0 ^ ^ ^ 1, so that the rest of the gravitational field 
is provided by the halo. For S12’s simulation, the fraction 
of active surface density was ^ = 0.5. The dependence of Q 
with radius is represented on figure 3. The scale invariance 
of this stability parameter leads as expected to a constant 
Q ^ 1.5 throughout most of the disc. It is only broken by the 
presence of the tapering functions which damp the surface 
density Et for the most inner and outer regions. 



J(l) 


Figure 2. Contours of the initial distribution function in action- 
space within the epicyclic approximation. The contours 

are spaced linearly between 95% and 5% of the distribution func¬ 
tion maximum. 


Q 



Figure 3. Variation of the local Q Toomre parameter with the 
angular momentum J^j^. It is scale invariant except in the in¬ 
ner/outer regions because of the presence of the tapering func¬ 
tions Tinner ^nd Touter - The unit systcm has been chosen so that 
Vo=G = Ri = l. 


3 THE SECULAR DIFFUSION EQUATION 

The figure 7 of S12 exhibits a ridge best seen in the con¬ 
tours of the distribution function in (Jr)—space. This ir¬ 
reversible diffusion feature of the DF was obtained through 
transient spiral features, while evolving a stationary Mes¬ 
tel disc sampled by pointwise particles. The formalism of 
secular diffusion should allow us to predict and explain its 
appearance. Let us first recall the main equations govern¬ 
ing secular diffusion. The secular dynamics intends to de¬ 
scribe the long-term aperiodic evolution of a self-gravitating 
system, during which small resonant and cumulative ef¬ 
fects can add up in a coherent way. These small effects, 
amplified through the self-gravity of the system, can be 
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seeded in two manners. As argued earlier, a possible first 
approach to describe such seeds is the secular diffusion equa¬ 
tion which describes the long-term evolution of a collision¬ 
less self-gravitating system undergoing external perturba¬ 
tions. The second approach is the Balescu-Lenard equation 
for inhomogeneous system which describes the evolution of 
closed self-gravitating system undergoing perturbations aris¬ 
ing from its own discreteness. In the first case, the pertur¬ 
bations are exterior and the system responds to it, whereas 
in the second case the perturbations are intrinsic and self- 
induced. Real discs respond to stimuli corresponding to a 
combination of these two formalisms. The finite number of 
stars in the disc, the presence of giant molecular clouds and 
of massive sub-halos around the disc all induce Poisson shot 
noise. Spiral arms in the gas distribution constitute another 
source of gravitational noise. Finally, the presence of a bar 
drives an additional coherent perturbation. The complex dy¬ 
namical history of a real stellar disc encompasses responses 
to all these stimuli. 

The context of the evolution of the Mestel disc stud¬ 
ied numerically by S12 corresponds to the Balescu-Lenard 
formalism. The perturbations originates both from the dis¬ 
crete sampling of the smooth DF, which is only represented 
by a finite number of particles, but also from numerical noise 
which can induce long-term and irreversible diffusion. The 
formalism of the Balescu-Lenard equation is much more in¬ 
volved than that of the external secular diffusion and will 
be the topic of future work. One can still however extract 
approximate yet interesting qualitative information for the 
long-term dynamics from the approach relying on external 
perturbation. In order to take into account the fact that the 
perturbations undergone by the system are created by the 
galactic disc itself, we will assume, as detallied later on, that 
the amount of noise a given location is generically propor¬ 
tional to the square-root of the local surface density, assum¬ 
ing that these intrinsic perturbations behave like a Poisson 
shot noise. 


3.1 The dressed secular equation in action space 


The distribution function introduced in equation (4) is a 
stationary distribution function, since it depends only on the 
actions coordinates J = (J(^, Jr)- The long-term evolution 
of the distribution in action-space is given by the secular 
diffusion equation, derived briefly in Appendix A (see also 
paper I for details), which takes the form 


dFo 

dt 




D, 




where the diffusion coefficients Dm{J) are given by 

Drr^iJ) = [[I-M]-'.e-[I-M]-'l . 
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(9) 


In this expression the response matrix, M, and the cross¬ 
power spectra C (see below) are functions of cj which should 
be evaluated at the resonant frequencies m fl, where the 
index m = {mr^nricf)) G corresponds to the Fourier coef¬ 
ficients associated to the Fourier transform with respect to 
the angles (Or, Off)) of the actions (Jr, Jcf))- To one specific m 
is associated one specific resonance. Throughout our calcu¬ 
lation, we will restrict ourselves to only three different reso¬ 


nances which are: the inner Lindblad resonance (ILR) cor¬ 
responding to m^^^) = (—1,2), the outer Lindblad 

resonance (OLR) for which = (1,2), and fi¬ 

nally the corotation resonance (COR) associated with cir¬ 
cular motion which satisfies = (0, 2). Indeed, 

S12 restricted disturbing forces to m(^ = 2, so that we may 
impose the same restriction on the considered azimuthal 
number Moreover, all the estimations presented in the 
upcoming sections have also been performed while consid¬ 
ering resonances with mr = ±2, which were checked to be 
subdmoninant, so as to justify our resonances restriction to 
the ILR, COR and OLR. Equation (9) for the diffusion coef¬ 
ficients also involves potential basis elements given by 
as introduced in the Kalnajs matrix method (Kalnajs 1976). 
Here (J) corresponds to the Fourier transform of index 
m with respect to the angles 9. Another key element of 
equation (9) is the response matrix M, which indicates how 
the system amplifies a given perturbation. It is given by 

, ( 10 ) 

m 

where one can see that the pole at the intrinsic frequency 
uj — ra Ft plays a crucial role for the amplification. 

In order to underline the physical meaning of these dif¬ 
fusion coefficients, one can rewrite them under the shortened 
form 

-^(a; = m.n), (11) 

|£(m,a;)| 

where qualitatively (see Chavanis (2012b) for the homoge¬ 
neous case), we have the following scalings 



where the coefficients b correspond to the basis decomposi¬ 
tion of the exterior perturbation so that 
These coefficients are then Fourier transformed with re¬ 
spect to time and one only needs to study their statis¬ 
tical ensemble-averaged cross-correlation defined in detail 
in equation (All). The diffusion coefficients are given by 
the ratio of the power spectrum of the external perturba¬ 
tions (|'0!^^(ct;) I ) divided by the gravitational susceptibility 
\s(m,uj)\‘^ of the disc. We will suppose that the exterior per¬ 
turbation, which represents the intrinsic noise of the system 
has a particular spectrum, since it originates from the galac¬ 
tic disc itself, so that 

(13) 

This assumption on the perturbations is relatively crude 
since we have only included a spatial dependence of the noise 
with J(^. For a system perturbed by a more realistic exte¬ 
rior source, the spectrum of perturbations is more coloured 
and depends on the full statistical properties of the exterior 
perturbers. Here the lack of dependence with the temporal 
frequency uo implies that the three resonances ILR, OLR 
and COR undergo the same perturbations at each position, 
even if they are not associated to the same frequencies of 
resonances m fl. As a consequence, there is no distinctions 
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between the perturbations felt at the inner and outer Lind- 
blad resonances. From the shape of the surface density in 
figure (1), one may see that the region of the inner taper¬ 
ing {J^ ~ 1.5) will be the most perturbed. With this as¬ 
sumption, expression (9) for the diffusion coefficients may 
be rewritten under the much simpler form 


Drr^iJ) 




(14) 


Equation (14) implies that the secular response of the system 
is the result of an arbitration between the system intrinsic 
noise and its local susceptibility. 

The diffusion equation (8) can be rewritten as the diver¬ 
gence of a flux in order to underline the exact conservation 
of total mass as 



We define as M{t) the mass contained in a volume V of the 
action-space at time t, so that 

M{t)= f d?JFo{J,t). (16) 

Jv 

Using the divergence theorem, the variation of mass due to 
secular diffusion can be seen as a flux of particles through the 
boundary S of this volume, with n being the corresponding 
exterior pointing normal vector, so that 

^ = f dS{m-n)DmiJ)m-^ . (17) 


In equation (17), the contribution from a given resonance m 
corresponds to a preferential diffusion in the direction m. 
This diffusion is anistropic in the sense that it is maximum 
for n oc ±m, and equal to 0, for n rn = 0. Two quantities 
influence the strength of the diffusion. On the one hand, 
the diffusion coefficient Dm{J) describes how sensitive the 
system is at a given location in action-space. On the other 
hand, the gradient m-dFoldJ describes how structured and 
inhomogeneous the system is at the same location. Such a 
formulation is of interest. First of all, it allows us to identify 
in each position Jr) what is the dominant resonance. 
It also permits us to identify the position of the peak of 
maximum flux, where the total flux, .7- tot, is defined as 



In this expression, the sum on the resonances m will be re¬ 
stricted to the ILR, OLR and COR resonances. From the 
contours maps of .T^tot, we are able to explain the ridge ob¬ 
served in S12 experiment. 

One may rewrite the diffusion flux from equation (18) 
using the slow and fast actions (Lynden-Bell 1979; 
Earn & Lynden-Bell 1996). We consider a given resonance 
m and introduce the change of coordinates 


F — 


J m 


jf _ J ^ 

^ nm. I I 


(19) 


where and are respectively the slow and 

fast actions associated to the resonance m = (mr, ?tI(^), 
— [nricf)^ —mr) corresponds to the direction perpendic¬ 
ular to the resonance, and |m| = ^/m rn. Using the chain 


A 



Figure 4. Variations of the response matrix eigenvalues A with 
the WKB-frequency kr- Red = small J^, Yellow = Large J^. 


rule, for any function X{J^, Jr), one has 

dX 


dX I I 

""a? ='"’I ajs 


( 20 ) 


We also naturally introduce the vector basis elements 
i^rn — = m'^/\m\) associated with this change 

of coordinates. In order to ease our qualitative description, 
we will assume that in the flux, equation (18), only one 
specific resonance m dominates. The diffusion flux J-rn, as¬ 
sociated to the resonance m, expressed within this rotated 
basis (J-^, Jfn), can now be rewritten under the form 

J^rr^iJ^,JL) = \mfDrr^iJ)■^etr,. (21) 

OJrn 


This rewriting shows that as soon as only one resonance 
dominates the secular dynamics, the diffusion flux will 
be aligned with this resonance, causing a narrow mono¬ 
dimensional diffusion in the preferential direction. Dur¬ 
ing this secular diffusion, particles conserve their fast action 
J^, which can henceforth be considered as adiabatically in¬ 
variant, whereas their slow action gets to change. We 
will show that such a mono-dimensional diffusion is indeed 
responsible for the ridge observed in S12 simulation. 

The evolution of a real disc is a combination of such 
resonant diffusions. Because stars are born on circular or¬ 
bits, action space is at first mostly populated close to the 
axis. Diffusion in the Jr—direction tends to increase the 
velocity dispersion within the disc and heats it. Diffusion in 
the J(^—direction brings stars from one nearly circular orbit 
to another of different radius and is called radial migration. 
This mechanism does not heat the disc, and because the 
density of stars does not change rapidly along the J^^—axis, 
it can go unnoticed. However, chemical evolution within the 
disc induce a radial gradient in metallicity with which radial 
migration can interact (Sellwood X Binney 2002). Indeed, 
near the Sun, it tends to wipe out the correlation between 
the ages and metallicites of stars, by bringing to the Sun 
both old metal-rich stars formed at smaller radii and young 
metal-poor stars formed at larger radii. 


3.2 The WKB tepid disc approximation 

One of the main difficulty of the implementation of the sec¬ 
ular diffusion equation (8) is to evaluate the diffusion co- 
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efficients, equation (9). This difficulty remains even within 
the assumption of a simple coloured noise, as introduced 
in equation (13). Indeed, it requires to define explicitly po¬ 
tential basis elements so that one can determine the 

response matrix from equation (10). The next step is to 
invert this response matrix M, in order to compute the 
diffusion coefficients Dm{J) from equation (9). To ease 
these calculations, one may rely on the WKB assump¬ 
tion (Liouville 1837; Toomre 1964; Kalnajs 1965; Lin & Shu 
1966), which assumes that the perturbations will take the 
form of tightly wound spirals, which in turn allows us to 
write Poisson’s equation locally. Such an approximation is 
well-suited to study S12’s experiment, because the secular 
evolution therein is sustained by the spontaneous genera¬ 
tion of transient spirals. Considering only such perturba¬ 
tions amounts to considering basis elements with specific 
properties and shapes. As explained in Appendix B, the 
main consequence of the WKB approximation is that the 
response matrix from equation (10) becomes diagonal. The 
expression of its eigenvalues, Xk^,kr,Ro{^)^ ^or a tepid disc 
reads 

_ _ C, ,\ _ 27vG^Ylt\kr \ ^ /'r)0^ 

Xk^,krRo{^) ~ X) 5 (22) 

where we have taken into account that only a fraction ^ of 
the disc is self-gravitating. Here kr corresponds to the radial 
frequency of the local spiral perturbation which is getting 
amplified, k^p is its azimuthal coefficient, which verifies kcp = 
2 for the ILR, OLR and COR, and Rq is the radius at which 
AC, St and x have to be evaluated. Here s is a dimensionless 
parameter reading 


LU k(f) 

s — - -—- . 

AC 

The dimensionless quantity x is given by 

CT kqr. 



(23) 


(24) 


Finally, the reduction factor, J^(s,x), (Kalnajs 1965; 
Lin & Shu 1966) is defined as 


-X T 

^ mr = l LFTJ 


4x] 


(25) 


Within the WKB approximation, one can show that the 
diffusion coefficients (fed by a stationary external perturba¬ 
tion -0^ oc depending only on position J(p, see equa¬ 

tion (13)) can be expressed in terms of the response eigen¬ 
values under the form 


Dm{J) — St 



dkr 



1 

1 — 


(26) 


where Jmr is the Bessel function of the first kind of index rrir 
and the integration on kr corresponds to an integration on 
all the radial frequencies of the tightly wound spirals, each 
one being amplified by the amplification factor 1/(1 —Afc^). 
Equation (26) is a novel result derived in Appendices A and 
C (see also paper I). The eigenvalues Xk^ have to be evalu¬ 
ated at the resonances so that uj = mfl. At such a resonance, 
one can see that for m = (mr^rricp), s is given by s = rrir. 
In order to handle the singularity of the equation (22) ap¬ 
pearing for s = ±l, one adds a small imaginary part to the 
frequency of evaluation, so that s^rrir + irj. Indeed, as long 


1 


1 '^max 



Figure 5. Dependence of the amplification factor 1/(1 — Amax) 
with the position for the ILR resonance (red) and the COR 
resonance (pink). Throughout the entire disc, the COR resonance 
is more amplified than the ILR and OLR resonances 


as r] in modulus is small compared to the imaginary part of 
the least damped mode of the disc, adding this complex part 
has a negligible contribution on the expression of Re(A). Fi¬ 
nally, we notice that for the Fourier coefficients associated 
to the inner/outer Lindblad resonances, one has = 
and 1. As equation (22) only depends on s^, these 

two resonances have the same response matrix eigenvalues, 
and as we have assumed that they are subject to the same 
noise, will therefore lead to diffusion coefficients of equal 
magnitude. 


3.3 Properties of the WKB equation 


Let us first study the behavior of the function kr i— Xk^ for 
a given resonance, m, and angular momentum, J(p. This 
function describes at a given position how much a pertur¬ 
bation with a frequency kr is locally amplified. Figure 4 
shows that for a given angular momentum there is a 
preferred frequency k^^^(J(p) for which X(kr) is maximum. 
One can now simplify the expression of the diffusion coef¬ 
ficients from equation (26) thanks to the behavior of the 
function kr ^ Xk^. Indeed, one can see that this function 
is peaked around k™^^{Jcp) with a typical spread equal to 
^kx{J(t))- Considering only the contribution from the region 
where Xk^ is maximum, equation (26) becomes 


Dm{J) — ^t(Rfir) A/ca sJ^rrir 
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2 


1-A 


max 


(27) 


The typical width of the amplification peak Akx is estimated 
via the width at half-maximum of the function kr i—A^^, so 
that Akx = ^ 1/2 ~ ^ 1 / 2 ? where ^ 1/2 fhe two so¬ 

lutions of the equation X{kr) = Amax/2. For the specific case 
of a Mestel disc considered by S12, the spread Akx and 
its position /cmax satisfy an additional property. Indeed, we 
have supposed that throughout the disc the radial-action 
spread was constant and we know from equation (3) 
that the epicyclic frequency k varies as l/Jcp. As a conse¬ 
quence, one has from equation (24) that x cx: {krJcp^. One 
may then rewrite the dependence with kr and J(p of the 
eigenvalues A under the form: X{kr,J(f)) = fiiJcp) f 2 {krJ(f)), 
where fi and /2 are given functions, which depend on the 
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resonance considered. Such a dependence of A with the radial 
frequency kr immediately implies that an additional scale- 
invariance property is satisfied so that Akx{J<i)) ^ 

/cmax oc XjJff. The inner regions of the disc have therefore 
larger eigenvalues spread than the outer regions. Such a de¬ 
pendence tends to enhance the susceptility of the most inner 
regions, which physically makes sense. 

Another important factor in the diffusion coefficients- 
given by equation (27) is the local amplification factor 
l/(l-Amax)^ which describes the strength of the ampli¬ 
fication due to the self-gravity of the system. As noted pre¬ 
viously, the ILR and OLR resonances possess the same am¬ 
plification factor. However, one can compare the strength of 
the amplification for the ILR and the COR resonance as seen 
on figure 5. One can note that the maximum amplification 
1/(1 —Amax) (~ 3 for the COR and ~ 1.5 for the ILR) re¬ 
main rather small so that the susceptibility of the disc is not 
too important. Moreover, one can note that the corotation 
resonance is more amplified than the inner/outer Lindblad 
resonances. However, what really represents the strength of 
the diffusion in the action diffusion map of the distribution 
is not the value of the diffusion coefficients Dm(«/) but the 
flux given by Dm(J) m\dF<^ldJ'). As the disc is tepid, one 
has for most of the regions \dF^ldJr \ ^ \dF^ldJcf\^ so that 
the ILR and OLR resonances for which vfir / 0 are favoured 
by the inhomegeneity of the distribution function compared 
to the COR resonances. This arbitration between the inho¬ 
mogeneity of the disc and its susceptibility determines the 
dominant regime of secular diffusion undergone by the disc. 


4 RESULTS 

In order to explain the ridge observed in S12, one can now 
compute on the plane (Jr, Jcf) the value of the diffusion co¬ 
efficients Dm{J) from equation (27) for the three main res¬ 
onances and obtain a numerical estimation of the maximum 
flux from equation (18), from which the secular diffusion will 
start. We will then vary the main features of the galaxy and 
its environment to trace its effects on the ridge. 

4.1 Reproducing the S12 experiment 

The results using the numerical prefactors from S12 are pre¬ 
sented on figure 6. The thin lines represents the contours of 
the diffused distribution function obtained by S12, which 
exhibits a narrow ridge of diffusion. Superimposed on these 
contours are shown the contours of the norm of the secular 
diffusion flux from equation (18), and the direction associ¬ 
ated to this flux. 

There is only one maximum peak of diffusion located 
in {Jr, Jcf) ^ (0.01, 1) from which the secular diffusion will 
unambiguously start. One should note that the predicted 
position of the peak of diffusion is slightly offset from the 
one observed in S12, which was around ^ 1.2. This dif¬ 
ference may be due to our crude model of intrinsic noise 
from equation (13), the fact that we are comparing the ini¬ 
tial diffusion direction to the final position of the ridge (see 
section 4.4 below), and/or possibly also to the limitations 
of the WKB approach which is only approximately accu¬ 
rate in a regime where the transients spirals are not suffi¬ 
ciently tightly wound. However, even so, the agreement on 


the regime of secular diffusion undergone by the disc remains 
quite good. To this maximum of the norm of the diffusion 
flux is also associated a direction of diffusion. The direc¬ 
tion of the ILR-diffusion is associated to the vector (—2,1) 
in the (J^)—plane, which makes an angle of 153° with 
the J(f)—axis. The direction of diffusion predicted within our 
approach is of approximately 120°. This quasi-alignement il¬ 
lustrates the fact that the ILR is the main resonance of the 
secular evolution of the tapered Mestel disc. S12 had found 
that the diffusion in action space was completely dominated 
by the ILR resonance, so that the ridge was aligned with the 
direction m = ttiilr. In our case, the slight misalignement 
observed has two origins. First of all, we considered a total 
secular flux, equation (18), summed on the three resonances 
ILR, OLR and COR, which all have different diffusion di¬ 
rections, so that it tends to slightly rotate the direction of 
the dominant resonance. 

Moreover, recall that our noise assumption from equa¬ 
tion (13) has no u = m-ft dependence, so that the ILR and 
OLR resonances possess the same susceptibility, which leads 
to an over-represent at ion the OLR resonance. However, we 
unambiguously recover that the S12 disc was in a ILR dom¬ 
inated regime, taking place in the inner regions of the disc. 
Since the ILR dominates the diffusion, let us use the analysis 
of mono-dimensional diffusion. The slow and fast actions as¬ 
sociated to the ILR resonance are given by J^lr oc J(/> — Jr/2 
and JjLj^ oc J^/2-\-Jr. In the neighbourhood of the diffusion 
peak, stars can therefore be assumed to diffuse along lines of 
constant along which their slow action J/lr changes. 

This leads to an mono-dimensional diffusion causing an heat¬ 
ing of the disc. 

Also, note from figure 6 that the diffusion flux norm is 
non-negligible only in a narrow band in J/lr- The size of this 
region will determine the typical width of the narrow ridge 
in the direction observed in S12 simulations. Starting 

from a narrow region in j(i^^ and diffusing predominantely 
in the J/lr— direction, one can therefore explain the ridge of 
limited extent in the (J^^, Jr)—plane observed numerically in 
S12, and correctly captured by the WKB limit of the secular 
diffusion formalism. 

This secular behavior of a typical Mestel disc dominated 
by the ILR resonance in the inner regions of the disc can be 
interpreted in the following way. Recall that the intensity 
of the secular diffusion is encoded by the total flux from 
equation (18). The use of a fraction, ^ = 0.5, for the surface 
density reduced the susceptibility of the disc and therefore 
reduced the amount by which perturbations can be ampli¬ 
fied through self-gravity. Consequently, the susceptibility- 
structure of the disc via Dm{J) is not the only key param¬ 
eter to determine the peaks of diffusion, but rather its in¬ 
homogeneous structure represented by the gradients of the 
distribution function dF^IdJ. This has two implications. 
Since the disc is tepid, the distribution function gradients 
are more important for the ILR and OLR resonances than 
for the COR, as resonances with a non-zero rur component 
are favored. Moreover, the gradients are the highest for the 
inner-regions, where the cut-off takes place, because of the 
presence of the tapering function Tinner (Jc^). This clarifies 
why it could have been expected that the secular diffusion 
would predominantly take place through a Jr—heating in 
the inner regions. 

An additional feature of such ILR-dominated diffusion 


8 Jean-Baptiste Fouvry, Christophe Pichon 



0 2 4 6 8 10 


Figure 6. Map of the norm of the total flux summed for the three resonances (ILR, COR, OLR) {hold lines). The contours are spaced 
linearly between 95% and 5% of the function maximum. The red vector gives the direction of the vector flux associated to the norm 
maximum (arbitrary length). The background contours correspond to the diffused distribution from S12 {thin lines)^ which exhibits a 
narrow ridge of diffusion. 



Figure 7. Dependence of the maximum peak of diffusion and its associated direction with the (ad hoc) properties of the perturbations of 
the system given by equation (28). The various curves follow flgure 6. From left to right: (rip, ap) = (1,1), (3, +oo), (3 , 1). Increasing 
the power index np tends to shift the peak position to higher whereas decreasing ap enhances the importance of the ILR and tends 
to align the direction of diffusion with miLR- 


is the typical temporal growth rate of the ridge in action 
space. As shown in Appendix C such an anisotropic and 
mono-dimensional diffusion, when started from a narrow hot 
spot, leads to a faster diffusion than the usual homogeneous 
heat equation. Indeed, the secular heating of a galactic disc 
through the ILR-resonance corresponds to a super-diffusive 
case, for which the diffusion of the distribution function dif¬ 


fuses does not follow the usual growth rate in \/t of the 
homogeneous heat equation: equation (C9) states that the 
scattering of the ridge in the direction grows linearly 

with the time t. 
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4.2 Modifying the properties of the perturbation 


The spectral properties of the noise play an important role 
in the detailled properties of the diffusion. In order to qual¬ 
itatively describe this impact, one can slightly modify the 
assumption (13) on the noise and consider the more general 
class of perturbation 






exp 


(Jp f^(t/^) 


(28) 


where rip is a power index introduced in order to modify the 
assumption on the Poisson shot noise, and ap is a dimension¬ 
less parameter (rescaled thanks to Q(J(^)) adding a depen¬ 
dence on the temporal frequencies u. Our initial assumption, 
equation (13) corresponds to the case rip = 1 and ap +oc. 
In order to illustrate the impact of the perturbation’s power 
spectrum on the properties of the secular diffusion, possi¬ 
ble ad hoc choices of the parameters (riprap) are shown on 
figure 7. Increasing the index rip enhances the relative im¬ 
portance of the surface density, and from the figure 1, one 
can see that it will favor the region around Jcf ~ 1.5, where 
the inner tapering takes place. The choices of the inner ta¬ 
pering function Tinner from equation (6) and its cutting scale 
Ri are then responsible for the location of the maximum sur¬ 
face density on figure 1 and therefore for the position of the 
peak of maximum diffusion. Moreover, adding a dependence 
on u breaks the degeneracy between the ILR and OLR res¬ 
onances. Indeed, the frequency associated to a resonance m 
is given by u = m ft. From expression (3) for the intrinsic 
frequencies, one can note that the frequencies of the three 
resonances satisfy the inequalities 


0 < CJILR < Ct;cOR < CJOLR • (29) 


Consequently, the addition of the decaying exponential in 
the noise model (28) tends to enhance the ILR resonance 
compared to the other resonances. In the sum in equa¬ 
tion (18) of the total flux, the vector contribution from the 
ILR dominates the other resonances, and we recover the fact 
that the direction of diffusion associated to the peak of secu¬ 
lar diffusion is closely aligned with ttiilr compared to what 
was observed on figure 6, with the simpler noise assumption 
given by equation (13), for which the COR resonance plays 
a more important role. Such ad-hoc experiments illustrate 
how the detailled properties and dependence of the pertur¬ 
bations with J(j) and uj can impact the characteristics of the 
secular diffusion. 

One could imagine situations where matching the ob¬ 
served secular response of families of disc to a given model 
for the external perturbation would provide means of char¬ 
acterising the statistical properties of their cosmic environ¬ 
ment. 


4.3 Modifying the tapering of the disc 

The inner tapering function Tinner from equation (6), which 
represents the bulge of the galaxy, plays a crucial role to 
secularly induce an ILR dominated peak of diffusion in the 
inner regions. Here Tinner is characterized by two parame¬ 
ters (z^, Ri), where v controls the sharpness of the tapering, 
whereas Ri is the scale at which it takes place. In order 
to illustrate the impact of Tinner on the secular diffusion, 
some modifications of Tinner are shown in figure 8. When v 


is increased, the surface density Et becomes steeper, so that 
the inhomogeneity of the system, represented by dTo/^«/, 
becomes more important. As a consequence, the peak of 
diffusion tends to migrate to the region of higher DF gradi¬ 
ents, which are in the vicinity of the scale radius Ri. On the 
other hand, decreasing v tends to enhance the importance of 
the susceptibility coefficients Dm{J) in the determination of 
the peak of diffusions. We noted in equation (27), that the 
scale-invariance properties of the Mestel disc impose that 
A/ca (xl/Jcf, so that the inner regions are naturally more 
susceptible. As a consequence, decreasing i/ tends to mi¬ 
grate the peak of diffusion inwards. Finally, as expected, 
modifying the scale radius Ri naturally leads to a similar 
displacement of the peak of diffusion. One should note that 
such modifications do not have any significant impact on the 
direction of diffusion. 


4.4 Modifying the temperature of the disc 

An additional way to modify the property of the disc is to 
change its temperature by varying the value of (7r, which 
encodes the radial-action spread of the distribution, as in 
equation (4). Its impact is illustrated in figure 9. Decreas¬ 
ing ar leads to colder discs, which tend to possess a wider 
peak of diffusion in regions slightly more external. A wider 
peak of diffusion will lead to a wider ridge when the secu¬ 
lar diffusion will take place. The impact of ar on the secular 
diffusion properties is in fact convolved, as it influences both 
the stability of the disc via Q, but also the detailed prop¬ 
erties of amplification of the disc, through the expression of 
the eigenvalues in equation (22) via y. Both the character¬ 
istics of the peak of diffusion and its direction of diffusion 
are therefore influenced by ar- 

Interestingly, as the ridge effectively increases locally 
the temperature of the disc (see Appendix C), Figure 9 sug¬ 
gests that it will in turn rotate the ridge in the ILR direction 
found by S12. Hence we can expect the discrepancy found 
on figure 11 to decrease in time through this process. 


4.5 Increasing the active fraction of the disc 

Let us now study one last feature of diffusion while modi¬ 
fying some of the characteristics of the disc and the halo. 
The behaviour of the function i—1/(1 —Amax) on fig¬ 
ure 5 showed that the COR resonance has higher ampli¬ 
fication factors than the ILR and OLR resonances. In order 
not to be dominated by the ILR resonance in the inner re¬ 
gions, one may try to increase the susceptibility of the disc 
so that the predominant diffusion will take place through 
the COR resonance. The eigenvalues of the response ma¬ 
trix from equation (22) can be increased via an increase of 
the active fraction of the disc. Figure 11 illustrates such 
changes. One can note in figure 11 that as ^ is increased, 
a significant COR-dominated region around the position 
(Jr, J<f) ~ (0 , 2) appears. This new diffusion region ends up 
being more important than the ILR-peak around the posi¬ 
tion (Jr, Jff)) = (0.01, 1), as illustrated on figure 10. Indeed, 
as ^ increases, both and A^ST^ increases, but since one 
has AmS<Am2^^<l, for A^a:^ close to 1, the COR reso¬ 
nance gets more amplified than the inner/outer Lindblad 
resonances. The fast and slow actions associated to the coro¬ 
tation resonance, are straightforwardly given by Jcor oc Jcf 
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Figure 8. Dependence of the maximum peak of diffusion and its associated direction with the inner tapering function Tinner from equa¬ 
tion (6). The various curves follow figure 6. Here Tinner is characterized by the pair (i/, Ri). The S12 case corresponds to (i/, Ri) = (4, 1). 
From left to right: (i/, = (3 , 1) , (5,1), (4,2). Reducing u reduces the sharpness of the tapering and therefore reduces the gradients 

of the distribution function so that the peak of diffusion migrates to the most inner regions. The scale radius Ri is as expected a crucial 
parameter to determine the position of the peak maximum. 




Figure 9. Dependence of the maximum peak of diffusion and its associated direction on the temperature of the disc, encoded by ar- The 
various curves follow figure 6. The S12 case corresponds to ar = 0.284. From left to right: ar = 0.20, 0.40, 0.50. Larger ar corresponds 
to hotter disc and therefore more stable. Colder discs have a diffusion peak with a wider extension and therefore will lead to wider 
ridges. 


and Jqq^oc Jr. As the fast action tends to be conserved 
during a secular diffusion dominated by only one resonance, 
we can conclude that the new peak of diffusion observed 
on figure 11 corresponds to the diffusion of nearly-circular 
orbits, which increase their angular momentum while 
their radial energy Jr remains small. As the active fraction 
^ is increased, we observe a transition from the ILR dom¬ 
inated heating of the inner region (J^^ ~ 1) to a regime of 
radial migration of quasi-circular orbits in more intermedi¬ 
ate regions (~ 2) as shown in figure 10. With such higher 
active fractions, the secular diffusion mechanism is now in a 
different regime of long-term evolution, mainly determined 
by the susceptibility of the disc via the diffusion coefficients 
rather than by the inhomogeneity of the distribu¬ 
tion function through dFe^jdJ. 


5 CONCLUSIONS 

The secular diffusion equation (Binney h Lacey 1988; 
Weinberg 1993; Pichon & Aubert 2006) of a self-gravitating 
collisionless system was re-derived and implemented in the 


WKB limit, using angle-action variables for tightly wound 
spirals in a tepid disc described by a Schwarschild distribu¬ 
tion function. In this limit, the functional form of the diffu¬ 
sion coefficient allowed us to identify the ridge found in ac¬ 
tion space by S12 for a stable Mestel disc. It originates from 
a resonant mono-dimensional diffusion, which is maximum 
at a specific locus in the inner regions of the disc. As the disc 
model investigated by S12 was somewhat singular, its global 
scale invariance is only broken by the addition of the inner 
and outer tapering functions (representing the bulge and the 
edge of the disc), which, as expected, also play an important 
role in the determination of the regime of secular evolution, 
hence the position of the peak of diffusion. The birth of a 
resonant ridge is therefore the result of a subtle fine tun¬ 
ing between many parameters of the system. Indeed, having 
for example a noise of the form (x 6 -d{uj — uJp), corre¬ 
sponding to a perturbation peaked at a specific frequency is 
not a mandatory condition to observe a resonant ridge. The 
self-gravity of the disc (via ^ and A), its susceptibility (via 
Dm{J)), its inhomogeneity (via dFo/dJ), its temperature 
(via (7^), its bulge structure (via Tinner), and the source of 
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Figure 10. For a given value of the active fraction ^ of the disc, 
one can identify a peak of flux associated to the ILR resonance, 
around (J^, J^) — (0.01, 1), and one to the COR resonance, close 
to (Jr, J(/>) — (0 , 2). This flgure represents the dependence of the 
norm of the peak flux from equation (18) for the two different re¬ 
gions of resonance, ILR (in blue) and COR (in purple), as the self¬ 
gravity of the disc is increased. For ^ ^ 0.68, the secular evolution 
of the disc becomes dominated by radial migration effects through 
the corotation resonance in the intermediate regions, rather than 
by an heating via the ILR resonance in the inner regions. 


perturbations (via all play a non negligible role in the 

creation and the properties of the resonant ridge, as shown 
in the ad-hoc experiments of the sections 4.2 to 4.5. 

The solar neighbourhood shows at least three indica¬ 
tions of the secular mechanism described in this paper: 
i) Groups of stars of a given age see their random ve¬ 
locities increase with the age of the group (Widen 1977; 
Aumer & Binney 2009). ii) Around the Sun, the velocity dis¬ 
tribution of stars is made of various streams of stars (Dehnen 
1998). Despite the fact that each stream is made of stars 
with different ages and chemistries, they respond to some 
perturbation in a similar way. iii) In the (J^^, Jr) space, one 
observes a depression of the density of stars near Jr = 0 and 
an enhancement for larger Jr, so that the disturbance in 
stellar density follows a curve consistent with resonant con¬ 
ditions (McMillan 2011). Given a detailed characterisation 
of the perturbations induced by e.g. the cosmic environment, 
one could study their effects on a typical self-gravitating col¬ 
lisionless galactic disc. In the context of the ongoing GAIA 
mission, this externally induced secular evolution is thought 
to be potentially a powerful approach to describe the long¬ 
term resonant radial migration of stars and its impact on 
the observed metallicity gradients (Sellwood & Binney 2002; 
Roskar et al. 2008; Schonrich & Binney 2009; Solway et al. 
2012; Minchev et al. 2013). 

More generally, the formalism of dressed secular diffu¬ 
sion fed by an exterior perturbations can be applied to any 
integrable model, and may lead to various secular diffusion 
scenarii depending on the structure of the galaxy and the 
properties of the spectrum of the external perturbations. 
The WKB formulation, when applicable, is very useful to 
yield a multiplicative amplification of the exterior pertur¬ 
bation and a tractable quadrature for the diffusion coeffi¬ 
cients, which under certain circumstances can be written 
algebraically (equation (27)). Such simplification provides 
useful insight into the physical processes at work, e.g. the 


relevant resonances, their loci and their relative strengths. 
Note that in principle, one could integrate the diffusion 
equation in time and show that secular evolution will drive 
the distribution function of the underlying disc into a state 
of marginal stability. Such iteration is postponed to another 
paper (Fouvry et al. 2014c, but see section 4.4 for a hint 
of the expected outcome). If we rid ourselves of the WKB 
tepid disc approximation, such equation could also possibly 
describe the secular diffusion of dark matter cusps in galactic 
centers via the external stochastic feedback processes within 
the inner baryonic disc. 

One of the limitations of the analysis presented in this 
paper is the simplified noise model used introduced in equa¬ 
tion (13), which does not depend on the temporal frequency 
uj nor the the radial frequency kr and is therefore only a func¬ 
tion of the position in the disc via Jcf. Such approximated 
perturbations intend to recreate the self-induced Poisson 
shot noise due to the discrete sampling of the smooth distri¬ 
bution function of the disc. A possible improvement would 
then be to extract the typical power spectrum of the per¬ 
turbations from numerical experiments, so as to use it as 
a refined seed for secular diffusion. In a companion paper 
(Fouvry et al. 2014a), we will explore the same WKB limit 
via the Balescu-Lenard equation. This perspective will en¬ 
able us to get rid of our crude noise approximation, since 
it will naturally encompass the self-indueed shot noise due 
to finite—iV effects and its impact on the secular diffusion 
of the quasi-stationnary distribution function, as long as we 
assume that transient tightly wound spirals are governing 
the evolution of the system. Such an approach will allow us 
to discuss quantitatively the expected timescales of secular 
evolution and the respective roles of the drift and diffusion 
components. In spite of this shortcoming, we have shown 
on this simplified experiment that the secular formalism de¬ 
scribes both the nature of the collisionless system (via its 
natural frequencies and susceptibility), and its nurture via 
the structure of the external power spectrum. Hence it pro¬ 
vides the natural framework in which to study the cosmic 
evolution of collisionless systems. 


Acknowledgements 

Many thanks to J. Binney, S. Prunet and P. H. Ghavanis for 
detailed comments. JBF thanks the GREAT program for 
travel funding and the department of theoretical physics, 
Oxford for hospitality. GP and JBF also thank the Institute 
of Astronomy, Gambridge, for hospitality while this inves¬ 
tigation was initiated. This work is partially supported by 
the Spin(e) grants ANR-13-BS05-0005 of the French Agence 
Nationale de la Recherche and by the ILP LABEX (under 
reference ANR-lO-LABX-63) which is funded by ANR-11- 
IDEX-0004-02. 


REFERENCES 

Aumer M., Binney J. J., 2009, MNRAS, 397, 1286 
Balescu R., 1960, Physics of Fluids, 3, 52 
Binney J., 2013, Dynamics of secular evolution, Falcon- 
Barroso J., Knapen J. H., eds., Gambridge University 
Press, p. 259 







12 Jean-Baptiste Fouvry, Christophe Piehon 



Figure 11. Map of the norm of the total flux summed for the three resonances (ILR, COR, OLR), when increasing the active fraction 
^ of the disc. The contours are spaced linearly between 95% and 5% of the function maximum for each case. From left to right: 
^ = 0.65, 0.68, 0.71. (Such values of ^ still comply with the constraint Q{J(p) > 1). 


Binney J., Lacey C., 1988, MNRAS, 230, 597 
Binney J., Tremaine S., 2008, Galactic Dynamics: (Second 
Edition), Princeton Series in Astrophysics. Princeton Uni¬ 
versity Press 

Born M., 1960, The Mechanics of the Atom. F. Ungar Pub. 
Co. 

Chavanis P.-H., 2012a, Physica A Statistical Mechanics and 
its Applications, 391, 3680 

Chavanis P. H., 2012b, European Physical Journal Plus, 
127, 19 

Dehnen W., 1998, AJ, 115, 2384 
Dubois Y. et ah, 2014, ArXiv e-prints 
Earn D. J. D., Lynden-Bell D., 1996, MNRAS, 278, 395 
Earn D. J. D., Sellwood J. A., 1995, ApJ, 451, 533 
Evans N. W., Read J. C. A., 1998, MNRAS, 300, 106 
Fouvry J. B., Piehon C., Chavanis P. H., 2014a, in prep 
Fouvry J. B., Piehon C., Prune! S., 2014b, submitted 
Fouvry J. B., et ah, 2014c, in prep 

Goldstein H., 1950, Classical mechanics. Addison-Wesley 
Heyvaerts J., 2010, MNRAS, 407, 355 
Jalali M. A., Hunter C., 2005, ApJ, 630, 804 
Jeans J. H., 1915, Monthly Notices of the Royal Astronom¬ 
ical Society, 76, 70 

Kalnajs A. J., 1965, Ph.D. thesis. Harvard University 
Kalnajs A. J., 1976, ApJ, 205, 745 
Kalnajs A. J., 1977, ApJ, 212, 637 
Lenard A., 1960, Annals of Physics, 10, 390 
Lin C. C., Shu F. H., 1966, Proceedings of the National 
Academy of Science, 55, 229 
Liouville J., 1837, ”j. math, pures appl.” 

Lynden-Bell D., 1979, MNRAS, 187, 101 
McMillan P. J., 2011, MNRAS, 418, 1565 
Mestel L., 1963, MNRAS, 126, 553 

Minchev L, Chiappini C., Martig M., 2013, A&A, 558, A9 
Nardini C., Gupta S., Ruffo S., Dauxois T., Bouchet F., 
2012, Journal of Statistical Mechanics: Theory and Ex¬ 
periment, 12, 10 

Piehon C., Aubert D., 2006, MNRAS, 368, 1657 
Piehon C., Cannon R. C., 1997, MNRAS, 291, 616 
Planck Collaboration et ah, 2013, ArXiv e-prints 
Polyanin A., 2001, Handbook of Linear Partial Differential 
Equations for Engineers and Scientists. Taylor & Francis 


Roskar R., Debattista V. P., Stinson G. S., Quinn T. R., 
Kaufmann T., Wadsley J., 2008, ApJ, 675, L65 
Schonrich R., Binney J., 2009, MNRAS, 399, 1145 
Sellwood J. A., 2012, ApJ, 751, 44 

Sellwood J. A., Athanassoula E., 1986, MNRAS, 221, 195 
Sellwood J. A., Binney J. J., 2002, MNRAS, 336, 785 
Solway M., Sellwood J. A., Schonrich R., 2012, MNRAS, 
422, 1363 

Toomre A., 1964, ApJ, 139, 1217 
Vauterin P., Dejonghe H., 1996, A&A, 313, 465 
Weinberg M. D., 1993, ApJ, 410, 543 
Weinberg M. D., 1998, MNRAS, 297, 101 
Weinberg M. D., 2001, MNRAS, 328, 321 
Wielen R., 1977, A&A, 60, 263 

Zang T. A., 1976, Ph.D. thesis. Massachusetts Institute of 
Technology 


APPENDIX A: SECULAR DIFFUSION 

Let us derive briefly the secular diffusion equation intro¬ 
duced in equation (8). For more details, see the companion 
paper Fouvry et al. (2014b). We consider a stationary dis¬ 
tribution function, Eo(J) (depending only on the actions J 
of the system (Jeans 1915)) undergoing an external pertur¬ 
bation. We also suppose that the gravitational potential of 
the system is fixed so that the action-coordinates mapping 
(x^v) ^ (0,«/), where 0 are the angles canonically associ¬ 
ated to the actions, does not depend on time. The distribu¬ 
tion function and the Hamiltonian of the system take the 
form 

F(J,0,t) = Fo(J,t) + /(J,0,t), 
iL( J, 6/, t) = iLo( J) + ^"(J, 6/, t) -h V^"( J, 6/, t), 

where / is the perturbation of the distribution function, 
is the perturbing exterior potential generated by an exter¬ 
nal system, and 'ijj^ is the self-response from the galactic 
disc generated via self-gravity. We assume that the pertur¬ 
bations are small so that f Fq and ^ '0o, where 

fjQ is the stationary background potential. We denote by 
ft = dHofdJ the intrinsic frequencies associated to the ac- 
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tions. Assuming that the disc evolves according to the colli¬ 
sionless Boltzmann equation, 


dt 


dF 

^ + («,F) = 0, 


(A2) 


with { , } a Poisson bracket, one obtains from equa¬ 
tion (Al), in the quasi-linear limit, two equations corre¬ 
sponding to the two timescales of the problem 

d{r+r) . 

dt^ ^ de dJ de 

< (A3) 

I dt (27r)2aj [J ^ do 


The first equation in (A3) describes the evolution of the 
perturbative distribution function / on the fast fluctuating 
timescale, whereas the second equation describes the long¬ 
term evolution of the stationary distribution function Fq 
in action-space. From the first equation, one obtains the 
expression of the diffusion coefficients appearing in equa¬ 
tion (8). We note by m = {mr^rricf) the Fourier coefficient 
associated to the Fourier transform with respect to the 
27r—periodic angles 0, so that the first equation of (A3) takes 
the form 

= 0. (A4) 

Following the matrix method (Kalnajs 1976), we introduce 
a biorthonormal basis of potential -0^^^ and densities 
satisfying 



Given this basis, the exterior potential -0^ and the self¬ 
potential may be written under the form 


< ■A ^ ^ (A6) 

^ (as). 

. P 

To shorten the notations, we also define the associated vec¬ 
tors a{t) = (ai(t), ..,ap(t),..) and h{t) = (5i(t),..., 5p(t),...). 
The next steps are to solve the equation (A4) for fm, then 
use the fact that the self-perturbing surface density ver¬ 
ifies that 5]^(cc) = f dv f{x,v) and recall that the trans¬ 
formation (cc,u) 1 -^ {0,J) is canonical so that it satisfies 
d^x d^v = d^O d^J. Given these remarks and assuming that 
dFo/dJ = cst. on the short timescale, one can show that the 
equation (A4) can be rewritten under the form 

a(ct;) = M(ct;) • |^a(ct;) + h(ct;)j , (A7) 

where the response matrix is given by equation (10). This 
expression describes how the perturbations get amplified on 
the short timescale. 

The next step is to use the second equation of (A3) 
to capture the secular evolution of the quasi-stationary dis¬ 
tribution function Fq in action-space. Introducing the sum 
of the two gravitational perturbations c(t) = a(t) + 6(t), one 
can show that the second equation of (A3) takes the form 


dFo 

dt 



m 


DmyJ-jtjTTl- 


(AS) 


where the diffusion coefficients are given by 


P,Q 


c’(i)rdre-“'”-“(‘-")cp(r). (A9) 

«/ —oo 


Using the matrix relation (A7), one can note that 
c = [I — M]“^' 6, so that equation (A9) can be expressed only 
as a function of the external perturbation b. The final step 
to derive the expression (9) of the diffusion coefficients is to 
introduce statistical averages, in order to consider only the 
mean response of a typical disc. We introduce the operation 
of ensemble average over many realizations of external per¬ 
turbations, noted as (.). When taking the ensemble average 
of equation (A8), one can assume that the response matrix 
M, the distribution function Fq and its gradient dFojdJ do 
not change significantly from one realization to another. In¬ 
deed, we intend to describe the effect of an averaged fluctua¬ 
tion on a given distribution function Fq representing a mean 
disc. Under these assumptions, equation (A8) becomes 


dFo 

dt 



m 


Dm(J,t) 


dFo 

' dJ 


(AIO) 


We finally introduce a stationarity hypothesis for the time 
evolution of the exterior perturbation on short timescales 
and therefore define the temporal auto-correlation of the 
exterior perturbation as 


Cki{ti -t2) = {hk{ti)lFi{t2)) . (All) 


From this definition, one can compute the ensemble average 
of the Fourier transformed term (^k to obtain 

(hk{io)h*{Lo')^ = 2 -k 5T>{LO-u2')Cki{>^), (A12) 


where C stands for the temporal Fourier transform of the 
autocorrelation function of the exterior potential. Injecting 
equation (A12) into equation (AIO) and using hermiticity 
arguments to show that only the real part of the diffusion 
coefficients matters, one obtains the final writing of the sec¬ 
ular diffusion equation, which reads 


dFo 

dt 




dFo 1 / 




(g)* 



[I-M]-'l (m-n) , 

J pq 


(A13) 


which is the same expression of the diffusion equation as 
introduced in equations (8) and (9). 


APPENDIX B: WKB COEFFICIENTS 


Let us derive briefly in this appendix the diffusion coef¬ 
ficients Dm{J) in the WKB limit. The key ingredient is 
to introduce specific basis elements well-suited to represent 
tightly wound spirals. We consider potential elements given 
by 




— 7 — Fail — 

(tTCT^) 


(R-Rof 


2cr2 


, (Bl) 


where the basis elements are indexed by three quantities. 
Here Ro is the central radius around which the Gaussian is 
centered, is an azimuthal number representing the angu¬ 
lar component of the basis elements, kr corresponds to the 
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radial frequency of the potential, and A = \/G/\kr\Ro is 
an amplitude tuned in order to ensure the correct normal¬ 
ization of the basis as imposed by equation (A5). Finally, 
cr is a scale-separation parameter ensuring the biorthog¬ 
onality of the basis elements. Introducing the typical ra¬ 
dial size of the system i^sys, one can show that under 
the assumptions of tight-winding krR^ ^ 1 and of scale- 
separation kr(T ^ Rsys/cr, the associated surface density ele¬ 
ments obtained via Poisson’s equation, are given 

by 

27rG 

In order to ensure the biorthogonality of the basis ele¬ 
ments, one can show that the distance between two basis 
elements and represented by /\Ro = RI — Rq and 
A/cr = kl — kr must satisfy 


ARo > cr or ARo = 0 , 
Akr ^ — or Akr = 0 . 


(B3) 


V cr 

With such explicit basis elements, one can compute their 
Fourier transforms with respect to the actions which read 




A 


*iTmr{,Rkr^ CXp 


(tTCT^) 1/4 

{R-Rof 


2cr2 


(B4) 


where J'mr is the Bessel function of the first kind of index 
TUr- Thanks to the WKB approximation which assumes that 
krRg ^ 1, the amplitude Hk^ and the phase shift 0% are 
given by 


Hkr ^ 


2Jr 


k'p 


Sr - -7r/2. 


(B5) 


Within this framework, one can now evaluate the response 
matrix elements from equation (10). Because of the assump¬ 
tions of tight winding, one can show that the response matrix 
is diagonal, so that two distinct WKB basis elements cannot 
interact one with another. Finally, we introduce the addi¬ 
tional assumption that the galactic disc considered is tepid, 
so that I dFo / d Jr \ ^ | dFo/dJ(f)\. Keeping only the dominant 
terms, the expression of the diagonal response matrix eigen¬ 
values for a tepid disc reads 


M 


27vG^^\kr\ , 

= —TTTZ -^ R{S, X) ■. 


[kl,k?,Ro],[kl,k‘‘,Ro] kl k? k2(1_s2) 


(B6) 


which is consistent with equation (22). Using the fact that 
response matrix is now diagonal, the expression (9) of the 
diffusion coefficients is easier to compute. 

The last step of the derivation is to express the ba¬ 
sis coefficients bp as a function of the perturbing exterior 
potential '0®^^ and to replace the sum on the basis index 
kr and Ro by continuous integrals, using Riemann formula 
^f{x)AxcB fdxf{x). As we have assumed that the ex¬ 
terior perturbation had the simple dependence from equa¬ 
tion (13), one finally obtains the expression of the diffusion 
coefficients in the WKB limit given in equation (26). 


APPENDIX C: GROWTH RATE OF RIDGE 

Let us estimate the rate at which the ridge in action space 
observed in S12 develops. As the evolution of the system 


is dominated by the ILR resonance, we introduce the as¬ 
sociated fast and slow actions from equation (19). In or¬ 
der to shorten the notations, we note the slow and fast ac¬ 
tions as J/lr = Js and — Jf- When considering only 
one resonance, the diffusion equation (15) becomes mono¬ 
dimensional and, for a given value of J/, up to constant 
prefactor in |m|, takes the form 


aup 

dt 


d 

dJs 


D{Js 


aup 

dJs 


(Cl) 


Equation (Cl) corresponds to a ID inhomogeneous heat 
equation. We suppose that initially the distribution function 
is concentrated within a narrow region in D, and we use a 
method based on self-similar solutions in order to estimate 
the typical growth rate of the ridge. Therefore, we intro¬ 
duce a self-similar Ansatz for the shape of the ID distribu¬ 
tion function of the form (see formula 1.3.6.8 from Polyanin 
(2001) and eq. (1.70) from Binney (2013)) 


F{Js,t) = YaP 



(C2) 


The coefficient a > 0 is for the moment unconstrained but 
will encode the speed with which the scatter of the distri¬ 
bution function increases. Indeed, one immediately obtains 
the time dependence of the standard deviation of the distri¬ 
bution function in the D—direction as 


o-Js = y {Ji )« p ■ (C 3 ) 


In order to be able to draw qualitative conclusions from a, we 
will assume that the anistropic diffusion coefficients satisfy 
a scaling property of the form 

D{k Js) = P D{Js), (C4) 


for any k > 0 and where 0 ^ 5 < 2 captures the structure 
of the anistropic diffusion. The homogeneous heat equation 
corresponds to the case D(Js) = cst, so that one has 5 = 0. 
Starting from the ansatz (C2) and using the rescaling change 
of variables K = JsjF one can show that the diffusion equa¬ 
tion (Cl) becomes 


a 

D+i 


T3c[A]+A 


dK 


,ab — 3 a ^ 


D{K) 


dFsc 

dK 


In order to have an equation valid for all time, both side 
must necessarily have the same time dependence, so that 
we obtain the condition 


ab — 2a—1. 


(C5) 


As a consequence, given the scaling coefficient 5, one con¬ 
strains a and then predicts the rate with which the standard 
deviation will increase during the diffusion. For example, 
for the homogeneous heat equation, one has 5 = 0, so that 
a = 1/2. Using equation (C3), one obtains that oc \/t, 
which is the usual growth rate of the scattering expected for 
an homogeneous diffusion. Moreover, one can also obtain 
the shape of the scale-invariant distribution function Dgcai 
through the ordinary differential equation 


d 


D{K) 


dFsc 

dK 


+ a 


Fsc[K] + K 


dFsc 

dK 


= 0 . 


(C6) 
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Using an ansatz of the form Fscai(^) = exp[—and as¬ 
suming that D(Js) = Do the solutions of the ordinary 
differential equation (C6) take the form 

K^-b 


Fsca,\{K) oc exp 


Do{2-hY 


(C7) 


In the case of the homogeneous diffusion 6 = 0, one ob¬ 
tains as expected a scale-invariant solution given by a Gaus¬ 
sian. We may now restrict ourselves to the case where 
the ILR resonace dominates, for which we want to esti¬ 
mate the temporal growth rate of an initially dense spot 
located in J^, on the Jr = 0 axis. Using the change of 
variables to the fast and slow actions, we introduce Jq and 
Jq such that (J^, Jr = 0) ^ (Jq, Jq )• As the diffusion takes 
place only in the J^—direction, we need to study the vari¬ 
ations of the diffusion coefficients given by the function 
js i^iLR(Jo+ Js 5 Jq)- We start from the expression of the 
diffusion coefficients obtained in equation (27). As the diffu¬ 
sion starts from Jr = 0, we may perform of limited develop¬ 
ment of the Bessel function at the origin, recalling that for 
X ^ 1, Jmr{x) oc SO that we obtain a dependence of 


the form 


Alr(Jo +isVo ) « j: 


(AfcA)^fc^ax 

(1 Amax)^ 




where the term between brackets has to be evaluated at the 
angular momentum Jcf corresponding to the fast and slow 
actions coordinates (Jq+Js, Jq)- This term tends a finite 
non-zero value for js ^ 0. As a consequence, for js ^ 0, one 
finally obtains the shortened dependence 

Drrv{Jo+j»,jf) OC Js . (C8) 


One can immediately conclude that for an ILR-dominated 
diffusion along its associated slow action, the scaling coeffi¬ 
cient of the anisotropic diffusion coefficients T>ilr is given by 
^iLR = 1. Using the relation (C5), one obtains that rilr = 1, 
so that the temporal growth rate of the standard deviation 
along the J/lj^— direction is given by 


oc t. (C9) 

This is the scaling presented in the main text. 








